BCR Work-Precision Diagrams
Samuel Isaacson and Chris Rackauckas
The following benchmark is of 1122 ODEs with 24388 terms that describe a stiff chemical reaction network modeling the BCR signaling network from Barua et al.. We use ReactionNetworkImporters to load the BioNetGen model files as a Catalyst model, and then use ModelingToolkit to convert the Catalyst network model to ODEs.
using DiffEqBase, OrdinaryDiffEq, Catalyst, ReactionNetworkImporters, Sundials, Plots, DiffEqDevTools, ODEInterface, ODEInterfaceDiffEq, LSODA, TimerOutputs, LinearAlgebra, ModelingToolkit gr() datadir = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr") const to = TimerOutput() tf = 100000.0 # generate ModelingToolkit ODEs @timeit to "Parse Network" prnbng = loadrxnetwork(BNGNetwork(), joinpath(datadir, "bcr.net")) rn = prnbng.rn @timeit to "Create ODESys" osys = convert(ODESystem, rn) u₀ = prnbng.u₀ p = prnbng.p tspan = (0.,tf) @timeit to "ODEProb No Jac" oprob = ODEProblem(osys, u₀, tspan, p) @timeit to "ODEProb DenseJac" densejacprob = ODEProblem(osys, u₀, tspan, p, jac=true)
Parsing parameters...done
Adding parameters...done
Parsing species...done
Adding species...done
Creating ModelingToolkit versions of species and parameters...done
Parsing and adding reactions...done
Parsing groups...done
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100000.0)
u0: 1122-element Vector{Float64}:
299717.8348854
47149.15480798
46979.01102231
290771.2428252
299980.7396749
300000.0
141.3151575495
0.1256496403614
0.4048783555301
140.8052338618
⋮
1.005585387399e-24
6.724953378237e-17
3.395560698281e-16
1.787990228838e-5
8.761844379939e-13
0.0002517949074779
0.0005539124513976
2.281251822741e-14
1.78232055967e-8
@timeit to "ODEProb SparseJac" sparsejacprob = ODEProblem(osys, u₀, tspan, p, jac=true, sparse=true) show(to)
──────────────────────────────────────────────────────────────────────────
──
Time Allocations
────────────────────── ─────────────────────
──
Tot / % measured: 329s / 99.0% 55.7GiB / 99.4%
Section ncalls time %tot avg alloc %tot a
vg
──────────────────────────────────────────────────────────────────────────
──
ODEProb DenseJac 1 257s 79.1% 257s 43.1GiB 78.0% 43.1G
iB
ODEProb No Jac 1 29.7s 9.12% 29.7s 5.20GiB 9.40% 5.20G
iB
ODEProb SparseJac 1 26.7s 8.20% 26.7s 5.26GiB 9.51% 5.26G
iB
Parse Network 1 8.08s 2.49% 8.08s 0.93GiB 1.68% 0.93G
iB
Create ODESys 1 3.46s 1.06% 3.46s 808MiB 1.43% 808M
iB
──────────────────────────────────────────────────────────────────────────
──
@show numspecies(rn) # Number of ODEs @show numreactions(rn) # Apprx. number of terms in the ODE @show numparams(rn) # Number of Parameters
numspecies(rn) = 1122 numreactions(rn) = 24388 numparams(rn) = 128 128
Time ODE derivative function compilation
As compiling the ODE derivative functions has in the past taken longer than running a simulation, we first force compilation by evaluating these functions one time.
u = copy(u₀) du = similar(u) @timeit to "ODERHS Eval1" oprob.f(du,u,p,0.) @timeit to "ODERHS Eval2" oprob.f(du,u,p,0.) # force compilation for dense and sparse problem rhs densejacprob.f(du,u,p,0.) sparsejacprob.f(du,u,p,0.) J = zeros(length(u),length(u)) @timeit to "DenseJac Eval1" densejacprob.f.jac(J,u,p,0.) @timeit to "DenseJac Eval2" densejacprob.f.jac(J,u,p,0.)
ERROR: syntax: expression too large
Js = similar(sparsejacprob.f.jac_prototype) @timeit to "SparseJac Eval1" sparsejacprob.f.jac(Js,u,p,0.) @timeit to "SparseJac Eval2" sparsejacprob.f.jac(Js,u,p,0.) show(to)
ERROR: syntax: expression too large
Picture of the solution
sol = solve(oprob, CVODE_BDF(), saveat=tf/1000., reltol=1e-5, abstol=1e-5) plot(sol, legend=false, fmt=:png)
